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ABSTRACT 

Using 3D global hydro simulations coupled with radiative transfer calcula¬ 
tions, we study the appearance of density waves induced by giant planets in 
direct imaging observations at near infrared wavelengths. We hnd that a 6Mj 
planet in a typical disk around a IMq star can produce prominent and detectable 
spiral arms both interior and exterior to its orbit. The inner arms have (1) two 
well separated arms in roughly m = 2 symmetry, (2) exhibit ~ 10 — 15° pitch 
angles, (3) ~ 180 — 270° extension in the azimuthal direction, and (4) ~ 150% 
surface brightness enhancement, all broadly consistent with observed spiral arms 
in the SAO 206462 and MWC 758 systems. The outer arms cannot explain obser¬ 
vations as they are too tightly wound given typical disk scale height. We confirm 
previous results that the outer density waves excited by a IMj planet exhibit low 
contrast in the IR and are practically not detectable. We also find that 3D effects 
of the waves are important. Compared to isothermal models, density waves in 
adiabatic disks exhibit weaker contrast in surface density but stronger contrast in 
scattered light images, due to a more pronounced vertical structure in the former 
caused by shock heating and maybe hydraulic jump effect. To drive observed 
pairs of arms with an external companion on a circular orbit, a massive planet, 
possibly a brown dwarf, is needed at around [r ~ 0.7", PA~ 10°] (position angle 
PA from north to east) in SAO 206462 and [r ~ 0.6", PA~ 10°] in MWC 758. 

Their existence may be confirmed by direct imaging planet searches. 
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disk interactions — radiative transfer — planets and satellites: formation 
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1. Introduction 


Detecting forming planets in protoplanetary disks is crucial for constraining the theory 
of planet formation, as it directly addresses two key questions: when and where do planets 
form. However, it is difficult to directly detect signals from forming planets in protoplanetary 
disks, and methods for indirect detection are necessitated. 

In protoplanetary disks, planets can excite spiral density waves through gravitational 
disk-planet interactions ( Kley fc Nelson]20T2 ). Density waves have much larger physical scale 
and can be much more prominent in observations than planets themselves. Once detected, 
they can be the smoking gun of the embedded unseen planets. Recently, high angular 
resolution direct imaging observations at near-infrared (NIR) wavelengths have found spiral- 


arm-like features in protoplanetary disks around young Herbig AeBe stars SAO 206462 (Muto 


et al.||2012 ; Garufi et al.||2013 ) and MWC 758 ( [Grady et al.||2013 ; Benisty et al.||2015 ). These 
observations obtained polarized intensity {7^, the linear polarization component 

in the scattered light) images as a way to suppress the central starlight (Hinkley et al.||2009 


Quanz et al.|2011 i |Hashimoto et al.||2011 ). These disks are relatively face-on, and show a pair 
of grand-design spiral arms often extending over 180° in the azimuthal direction. These arms 
are located at several tens to a couple of hundred AU from the center, and they are nearly 
in m = 2 rotational symmetry. At this moment resolved mm observations of these objects 
suffer from insufficient spatial resolution, and could not pinpoint the mm counterparts of 


these arms (Ghapillon et al. 2008; Isella et al. 2010; Andrews et al. 2011; Isella et al. 2013 


Perez et al. 2014; Marino et al. 2015). Also, the parent disks are recognized as transitional 


disks (Espaillat et al. 2014), with a giant cavity at the center. 


The origin of the observed spiral arms is unclear, though they have been widely specu¬ 
lated to be density waves excited by embedded unseen planets. In linear theory of density 
waves, the pitch angle, defined as the angle between the azimuthal direction and the tan¬ 
gent of the waves, is set by the aspect ratio h/r of the disk, thereby the disk temperature 
(Rafikov 2002; Muto et al. 2012). Fitting the shape of observed arms to linear theory and 


assuming that the planets are interior to arms has sometimes led to unphysically high disk 
temperature ( Benisty et al.||2015 ). Recently, in a pioneering study Juhasz et al. (2015) calcu¬ 
lated the surface density structures of waves using 2D locally isothermal hydro simulations. 
The resulting disk structures were puffed up in the vertical dimension assuming a Gaussian 
profile, and synthesized NIR images were produced using 3D radiative transfer simulations. 
It was concluded that the amplitude of the waves induced by a IMj planet at 25 AU in a 
typical disk is too small to be visible with current NIR facilities. 


In this paper following Dong et al. (2014), we revisit the observational appearance of 


planet-induced density waves in NIR imaging observations. The primary question we address 
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is whether planets can drive density waves that resemble the directly imaged spiral arms. 
To do this, we carry out 3 D global hydrodynamics simulations using the Athena++ code to 
calculate the density structures of disks. Both locally isothermal (ISO) and adiabatic (ADI) 
equations of states are explored, and we highlight both the inner and outer density waves. 
The resulting 3 D hydro models are fed into 3 D Monte Carlo Radiative Transfer (MCRT) 
simulations to produce synthetic model images at 1.6/im {H band) using the ([Whitney et al. 


2013) code. The hydro and MCRT simulations are introduced in Section the modeling 


results are presented in Section followed by a short summary in Section 


2. Simulation Setup 


2.1. Hydro dynamical Simulations 


The detailed discussion of the hydro simulations using the Athena++ code and the den¬ 
sity wave structure will be presented in Zhu et al. (in prep). Here we provide a short summary 
of the models. Athena++ is the newly developed grid based code using a higher-order Go¬ 
dunov scheme for MHD and the constrained transport (CT) to conserve the divergence-free 


property for magnetic fields (Stone in prep.). Compared with its predecessor Athena (Gar- 
diner fc Stone|20d5 i 2008; Stone et al.|2008[ ), Athena-|--|- is highly optimized for efficiency and 
uses flexible grid structures, allowing global numerical simulations spanning a large range 
of scales. We have carried out global 3-D simulations in spherical coordinates r (radial), 0 
(polar), and (f) (azimuthal) with both isothermal and adiabatic equation of states. In the 
adiabatic runs the adiabatic index 7 is chosen to be 1.4. 


The initial disk temperature T is constant on cylinders 

-1/2 

r(r,z) = To(-' 

Vro 


( 1 ) 


Therefore the initial disk scale height h = Cg/D oc 


oc r 


1.25 


(cg is the sound speed 


and D is the angular frequency in the disk), and h/r = 0.1 at 50 AU. The disk density and 
azimuthal velocity are set to maintain hydrostatic equilibrium (e.g. Nelson et al.||2013 ). The 
initial disk surface density Sq decreases as r~^ in all our simulations. 


The grids are uniformly spaced in log r, 9, (p with 256 x 128x688 grid cells in the domain 
of [15, 150 (AU)] X [7r/2-0.6, 7r/2-|-0.6 ]x[0, 27r]. Gonstant a viscosity with a = 10“^ has been 
applied. At inner and outer boundaries, all quantities are fixed at the initial states. We have 
run the simulations for 20 planetary orbits. The planet is on a flxed circular orbit, and there 
is no accretion from the disk onto the planet. The time scale is chosen so that density waves 
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have established throughout the whole disk (longer than the sound crossing time) while a 
big and deep gap has not been opened (and consequently no vortex generated at the gap 
edge). The gap opening process depends on the level of turbulent stress in the disk, which 
is poorly constrained. In total, we carry out 6 models, and their properties are summarized 
in Table [H 


2.2. Monte Carlo Radiative Transfer Simulations 


Density distributions obtained in our hydro simulations were post-processed via the 


3D MGRT calculations using the code dev 
used to model protoplanetary disks in tl 

eloped by Whitney et al. 

to 

o 

CO 

), which has been 

le past (e.g. 

Hashimoto et al. 

2012 ; 

Zhu et al. 

2012 ; 

Dong et al. 


Follette et al. 

2013; Grady et al.||2013; |Dong et al. 

2014 

). In 


MCRT simulations, photons from the central star are absorbed/reemitted or scattered by 
the dust in the surrounding disk. The temperature in each grid cell is calculated based 


on the radiative equilibrium algorithm described in Lucy (1999). The anisotropic scattering 


phase function is approximated using the Henyey-Greenstein function (Henyey & Greenstein 


1941). Polarization is calculated assuming a Rayleigh-like phase function for the linear 


polarization (White 1979). Full resolution synthesized PI images at H band (1.6/im) are 
produced for all models with a pixel size of 0.7 AU/pixeQ These images are then convolved 
by a Gaussian point spread function (PSF) with a full width half maximum (FWHM) of 0.06" 
and assuming a distance of 140 pc, to achieve an angular resolution comparable with NIR 
direct imaging observations using Subaru, VLT, and Gemini (the FWHM of a theoretical airy 
disk is 1.028A/D ~ 0.04" for a primary mirror with a diameter D = 8.2matA = 1.6 /rm). 


The MGRT simulation setup is similar to Dong et al. (2014). We construct a 3D disk 


structure in spherical coordinates with the same grid structure as in the hydro models (0 = 0 
is west in all images). The central source is a typical Herbig Ae/Be star with a temperature 
of 10^ K and a radius of 2Rq. All simulations are run with 3 billion photon packages. The 
dust grains in the disk are assumed to be interstellar medium (ISM) grains ( Kim et ^[1994 ) 
made of silicate, graphite, and amorphous carbon. Their size distribution is a smooth power 
law in the range of 0.002-0.25 pm followed by an exponential cut off beyond 0.25 pm. These 
grains are small enough so that they are dynamically well coupled to the gas, and thus have a 
volume density linearly proportional to the gas density. The optical properties of the grains 


can be found in Figure 2 in Dong et al. (2012b). The total mass of the ISM grains is assumed 


^In this work, the physical quantity recorded in all model images is the specific intensity in unit of 
[mJy arcsec“^], or [ergs s“^ cm“^ Hz“^ arcsec”^]. 
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to be 2 X 10-^M(£1 


3. Modeling Results 

In this section we present the modeling results. The artificial central cavity in the 
hydro models produces a bright rim structure at its edge in scattered fight images, which 
contaminates disk structure within ~ 35 AU from the center in convolved images. Therefore, 
we focus only on the outer arms in the fp = 50 AU models and the inner arms in the 
Tp = 125 AU models, which are all located at r > 50 AU. 


3.1. Density Structure of the Models 


The surface density perturbation S/Sq (where Eg is the initial surface density) of the 
hydro models are shown in the left column in Figures and When Mp/M* ^ (h/r)^, 
a planet excites one inner and one outer waves originated from the planet (the primary 
arms from now on) and another two shifted by roughly 180° in the azimuthal direction 
(the secondary arms from now on)j^ while their strengths depend on Mp, /j/r, and EOS. 
The primary waves are always stronger than the secondary waves, which are practically not 
recognizable in 1ISO50 and 1ADI50 runs but visible in all other cases. In IMj models, S^ave 


never exceeds Eq by 50%, consistent with Juhasz et ah (2015), while in the 6Mj models 


the surface density enhancement of the primary waves can be more than 100% (Table [^. 
E enhancement in ISO models is close to twice that in the corresponding ADI models. We 
note that since in 3D the disk is not hydrostatic in the vicinity of the shock, its vertical 
structure cannot be fully described by puffing up the surface density distribution from 2D 


simulations assuming hydrostatic (Gaussian) vertical density profile, as performed in Juhasz 


et ah (2015). 


^This corresponds to, for example, a total gas disk mass of 0.02 Mq, a 100:1 gas-to-dust-mass-ratio, a 
10% dust mass fraction in the small ISM grains, and the rest 90% in the large grains that have already 
settled to the disk midplane and do not affect NIR scattering. 


^The two spiral arms are seen in previous simulations, e.g. Figure 10 of 
commented in|Juhasz et al.|(2015). 


de Val-Borro et al. 


(2006), and 
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3.2. Scattered Light Images of the Outer Arms 


Full resolution and convolved H band images for all models are shown in Figure and 
The outer arms in the IMj models are practically not traceable in both full resolution images 
and convolved images, echoing Juhasz et ah (2015) (they are recognizable in the 1/r^-scaled 


convolved images as the dynamical range of the background is largely suppressed). The peak 
surface brightness enhancement along the waves is r\j 10 — 20% (Table [^, producing only 
marginal flnctnations on radial and azimuthal image prohles (Figure]^. The 6 Mj models 
produce much more prominent arms (Figure and |^. The primary arm in 6ADI50 can be 
up to ~2 times brighter than the background, and it is radially extended, with a radial width 
~ 0.2" at r = 0.5". Note that the contrast of arms in our convolved images are expected to be 
higher than in observations due to (1) model images are convolved by a Ganssian PSF, which 
has a similar kernel but less extended wings than a realistic PSF, and ( 2 ) no observational 
noises and instrument effects are included, which will further weaken the contrast of features. 
The 2 outer arms have a rough m = 2 symmetry, and are tightly wound, with pitch angles 
no bigger than 6 ° (Table , which is roughly half the values in observed systems (listed in 
Table as well for comparison). They may not be distingnishable from a ring structure in 
observations with low S/N. 

An interesting finding is that arms in ADI models are more prominent than in ISO 
models. The arm snrface brightness enhancement in 1ISO50 is only half of that in 1ADI50, 
despite having a snrface density enhancement 2 times larger. 6 Mj models show similar 
effects. This is caused by a combination of shock heating and hydraulic effects in waves. As 


predicted by 

Goodman & Rafikov 

( 200 : 

U; 

Rafikov 

( 200 ^ 

1 ) and shown 

in simulations ( 

Dong 

et al. 

2011 ; 

Dnffell & MacFadyen ! 

2013; 

Zhu et al. 

2013; 

Richert et al. 

2015), density 

waves 


shock at a conple of scale heights away from the planet’s orbit when Mp > M*(/i/r)^. Shock 
heating heats material along the waves in ADI cases, resulting in a higher scale height than in 
ISO cases and puffing up the waves in the vertical direction {h = Cg/O oc \/T). The hydraulic 
effects seen in 3D simulations of wave shocks by Boley & Durisen ( |2006| ) also makes ADI 
waves more pnffed np. Also interestingly, in the Mp = 6Mj models while both arms are 
visible in 6ISO50 with roughly equal strengths (also note the brightening behind the planet, 
Jang-Condell & Turner 2012 ), in 6ADI50 the primary arm is much more prominent. This 


may indicate a difference in the excitation and dynamics between the two arms, as the shock 
heating and hydraulic effects are more effective in the primary arm. 

Lastly, we perform experiments in which we collapse our 3 D hydro disks to 2 D, then 
puff them up back to 3D assuming Ganssian profiles in the vertical structure, similar to 


what was done in Juhasz et al. (2015). Results show that a full 3D hydro treatment does 


not affect the detectability in the outer arms excited by a IMj planet, however in Mp = 6 Mj 
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models it does produce noticeable dilferences in images, particularly for the inner arms. A 
more detailed discussion can be found in Zhu et al. (in prep). 


3.3. Scattered Light Images of the Inner Arms 


In our simulations the inner arms are open and prominent. In convolved images, the 
arms are up to 30 — 40% brighter than the background region at the same radius in 1IS0125, 
and 150 — 160% brighter in 6IS0125 (Table [^. Interestingly, although the secondary arm 
are much weaker than the primary in surface density contrast, the two are roughly equally 
bright in scattered light images. In the azimuthal direction, the primary and secondary arms 
extends over ~ 180° and ~ 270°, respectively. The 2 arms are roughly 180° rotationally 
symmetric (in both location and contrast). The pitch angle is about 15° for the primary arm 
and 10° for the secondary. The typical FWHM of the arms in the radial direction is ~ 0.1". 


In general, the morphology of the inner arms is quite similar to the observed ones 
in SAO 206462 (Muto et al. 2012 | Garufi et al. 2013) and MWC 758 ( |Grady et al. |201^ 


Benisty et al. 2015), in terms of (1) the rough m = 2 symmetry, (2) pitch angle, and (3) 


azimuthal extension. We note that the shape (i.e. pitch angle) of these inner arms cannot be 
described by the linear theory in Rafikov (2002); Muto et ^ (2012), as the waves are highly 
nonlinear. A comparison between MWG758 and a rescaled 6IS0125 model with the actual 
source distance and angular resolution is shown in Figure]^ In addition, for observed arms 
while it is difficult to measure the surface brightness enhancement (difficulty in defining a 
true “background”) and radial width (noises in the data) in an equally clean way as we do to 
our theoretical models here, observed arms are mostly in the bulk part 100 — 300% brighter 
than the background and around 0.1" in width in the radial direction, broadly consistent 
with the 6IS0125 model. In general, the morphology of the inner arms in the 6IS0125 run 
is tantalizingly similar to the observed ones. 


4. Summary 


Our main conclusions are: 


1. The inner spiral arms excited by a massive planet {q = Mp/M* = 6 x 10“^) can 
be visible in scattered light images under current NIR direct imaging capability, and 
they have a morphology resembling the observed spirals, such as in the SAO 206462 


(Muto et al. 

2012 [ |Garufi et aL| 

2013 

) and MWC 758 ( 

Grady et al. 

2013[ |Benisty 


et al. 2015) systems (Figure]^. They exhibit (1) roughly m = 2 symmetry, (2) pitch 






























angles in between 10 — 16°, (3) an azimuthal extent of ~ 180 — 270°, and (4) ~ 150% 
surface brightness enhancement relative to the background, all broadly consistent with 
observations. Note that the shape of these arms cannot be fully described by the 


weakly non-linear theory (Rafikov 2002). 


2 . The outer spiral arm excited by a g = 6 x 10“^ planet can also be visible. However they 
cannot explain the observed spiral arms, as they are too tightly wound given typical 
disk scale height (pitch angle < 4°, 2.5 times or more smaller than observations), 
which has been noted in the case of MWC 758 by Benisty et al. ( 2015| ). In addition, 
the contrast of the isothermal waves is lower than observed spiral arms, while the more 
prominent ADI waves are too extended in the radial direction, with a radial width a 
factor of 2 or more larger than observations. 


3. Planet induced density waves have a higher surface density enhancement in isothermal 
disks than in adiabatic disks. However, in scattered light images the arms in adiabatic 
disks are more prominent, thanks to a more extended structure in the vertical direction 
due to shock heating and possibly hydraulic effects. Also, 3D effects in isothermal waves 
can be important (Zhu et al. in perp.) 


One key factor in the appearance of planet-induced density waves may be that the gap 
cannot be fully opened by the planet given the duration of our simulations, so that the bulk 
part of waves are not located in a deep gap region. In our models, the planet is located 
roughly at a distance 1.5-3 times the main part of the inner arms. If observed arms are 
driven by a companion on a circular orbit outside the arms, then a ~ lOMj companion 
(possibly a brown dwarf) is needed at around 100 AU in SAG 206462 (at a distance of 
140 pc, |Muto et al.||2012| ), and 160 AU in MWC 758 (at a distance of 280 pc, Grady et al. 


2013). They can be excellent targets for direct imaging observations.. However, the waves 
may also be excited by a recent ffyby, or a companion on a very eccentric orbit. Lastly, 
it has been speculated that the cavities in transitional disks are opened by multiple giant 
planets. Observed global-scale m = 2 arms so far are all found outside the central cavities 
in transitional disks. The hypothetical outer companion suggested in our work cannot be 
associated with the putative cavity-opening planets; the latter would reside inside the cavity 
and their density waves are unlikely to be detectable. 
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Table 1: Model Properties 


Model/Gbservation 

EOS 

Q 

I'p 

(^wave,peak ^o)/^0 

(-Iwave,peak -Ib)/-Ib 

Radial FWHM 

Pitch Angle 



xlO-3 

AU 

% 

% 

arcsec 

degree 

(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

1ISO50 

ISO 

1 

50 

49(OP) 

lO(OP) 

N/A 

6(OP) 

1ADI50 

ADI 

1 

50 

28(OP) 

20(OP) 

N/A 

6(OP) 

1IS0125 

ISO 

1 

125 

31(IP), 21(IS) 

40(IP), 30(IS) 

~ 0.09(IP), ~ 0.08(IS) 

16(IP),10(IS) 

6ISO50 

ISO 

6 

50 

150(OP),61(OS) 

N/A 

~ 0.09(OP), ~ O.l(OS) 

4(OP),2(OS) 

6ADI50 

ADI 

6 

50 

77(OP), 36(OS) 

190(OP) 

- 0.24(OP) 

4(OP) 

6IS0125 

ISO 

6 

125 

125(IP), 48(IS) 

160(IP), 150(IS) 

~ O.l(IP), - 0.08(IS) 

15(IP),10(IS) 

SAO 206462*^ 





^ 100 (east), ^ 300 (west) 

~ 0.11 (east), ^ 0.13(west) 

11 (east), 10 (west) 

MWC 758^ 





^ 200 (east), ^ 250 (west) 

^ 0.07(east), ^ 0.1 (west) 

11 (east), 10 (west) 


Note. — Properties of the models. (1) Model name. (2) Equation of state in the hydro simulations: 
isothermal (ISO) or adiabatic (ADI). (3) q = For systems around a IMq star the list number is 

just planet mass in Mj. (4) We only look at outer arms in the = 50 AU models and inner arms in the 
Tp = 125 AU models. (5) Peak surface density enhancement on the waves. The letters in parenthesis are for: 
(O) - outer arms; (I) - inner arms; (P) primary arms (originated from the planet); (S) - secondary arms. 
The secondary arm in 1ISO50 and 1ADI50 are not recognizable. (6) Peak surface brightness enhancement 
on the waves in convolved images. Jb is the surface brightness in the background region around position 
angle b = 270 — 325° at the same distance from the center as the arms (we avoid the dark shadowed regions 
right behind the arms). The two outer arms in 6ISO50 are both prominent and tightly winded up, which 
makes it difficult to define a “background” at the distance of the arms. The secondary arm in 6ADI50 
cannot be easily distinguished from the background. (7) FWHM of the main part of the arms in the radial 


to 

I 


direction. (8) Averaged pitch angle of the arms, (a) The more recent and sharper VLT images in Garufi 


et al. ( 2013j SAG 206462) and Benisty et al. ( |2015 MWC 758) are used for the pitch angle and radial 
width measurements, while older Subaru datasets in Muto et al. | (2012, SAG 206462) and Grady et al. (2013 


MWG 758) are used for surface brightness enhancement measurements as the absolute polarized intensity 
calibration is not always available in the new VLT datasets ( Benisty et al.||2015) . We only focus on the 
2 major arms on the east and west sides in each case. It is sometimes difficult to perform measurements 
on observed arms in an equally clean way as we do to our models due to observational noises and time 
variations. 
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S/So band PI/(mJy arcsec iJ band PI/(mJy arcsec 

0 0.5 1 1.5 2 1 3 10 30 0.5 1 1.5 



-100 0 100 -100 0 100 -100 0 100 


X (AU) 


Fig. 1.— Model results of the Mp = IMj models, showing the relative surface density 
perturbation S/Eq (left column, linear scale); H band (1.6 fim) fnll resolution polarized 
intensity images (middle, logarithmic scale, the central 0.2" has been masked out to high¬ 
light the structures in the outer disk); and 1/r^-scaled convolved images (right, linear scale, 
convolved by a Ganssian PSF with FWHM= 0.06"). Note that gaps have not been fully 
opened by the planets. 
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Fig. 2.— The same as Figure but for the 6 Mj models. 
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Fig. 3.— Top: Radial profiles of convolved images at 90°, 180°, and 270° away from the 
planet, as illustrated in the cartoon. Bottom: Azimuthal profiles of convolved images at 
r = 0.5" (70 AU) from the star. The (dots, triangles) mark the intersections with the 
(primary, secondary) arms. 














16 



Fig. 4.— Comparison between 6IS0125 (left, 1/r^-scaled convolved image) and MWC 758 
(right, Benisty et al. 2015). The model is rescaled to fp = 158 AU, and the Gaussian PSF 
convolution adopts the source distance (280 pc) and the angular resolution of the observation 
(FWHM=0.03"). Units are arbitrary. The green dot in the model image marks the location 
of the planet (Mp/M* = 6 x 10“^). 











